In [54]:
%matplotlib inline
import datacube
from datacube_stats.statistics import GeoMedian
import sys
sys.path.append('../Scripts')
from dea_datahandling import load_ard
from dea_plotting import display_map
from dea_plotting import rgb
In [55]:
dc = datacube.Datacube(app='Geomedian_composites')
In [56]:
# Set the central latitude and longitude
central_lat = -28.200
central_lon = 153.169
# Set the buffer to load around the central coordinates - max 0.1 deg
buffer = 0.1
# Compute the bounding box for the study area
study_area_lat = (central_lat - buffer, central_lat + buffer)
study_area_lon = (central_lon - buffer, central_lon + buffer)
In [57]:
display_map(x=study_area_lon, y=study_area_lat)
Out[57]:
In [58]:
# Create a query object
query = {
'x': (central_lon - buffer, central_lon + buffer),
'y': (central_lat - buffer, central_lat + buffer),
'time': ('2019-10-01', '2019-12-31'),
'measurements': ['nbart_green', 'nbart_red', 'nbart_blue'],
'output_crs': 'EPSG:32756',
'resolution': (-10, 10),
'group_by': 'solar_day'
}
# Load available data
ds = load_ard(dc=dc,
products=['s2a_ard_granule', 's2b_ard_granule'],
**query)
# Print output data
print(ds)
In [59]:
# Set the timesteps to visualise, this
timesteps = [1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23]
# Generate RGB plots at each timestep
rgb(ds, index=timesteps)
In [60]:
# Compute geomedian using all observations in the dataset
geomedian = GeoMedian().compute(ds)
In [61]:
print(geomedian)
In [62]:
# Plot the result
rgb(geomedian, size=10)
